#################################################################
## Purpose: REPLICATION for "Criminal Victimization and Agency Attitudes in Mexico"
## Author: Cassy Dorff, 2024
## Step 2 cont'd: Tables for all models
#################################################################

#################################################################
rm(list=ls())
source(paste0("~/Desktop/replication/0_setup.R"))

# load the data
load(file=paste0(pathData,"1_nationalsample_data.rda"))
load(file=paste0(pathData,"2_oversample_data.rda"))

# load models data saved from 2_models_run.R
load(file=paste0(pathData,"model1_results.rda"))
load(file=paste0(pathData,"model2_results.rda"))
load(file=paste0(pathData,"model3_results.rda"))
load(file=paste0(pathData,"model4_results.rda"))
load(file=paste0(pathData,"model5_results.rda"))
load(file=paste0(pathData,"model6_results.rda"))
load(file=paste0(pathData,"model7_results.rda"))
load(file=paste0(pathData,"model8_results.rda"))
load(file=paste0(pathData,"model9_results.rda"))
load(file=paste0(pathData,"model10_results.rda"))

# short names
ndata <- national
#################################################################

#################################################################
# DATA SET UP
# set ref to largest group (vote)
ndata$ref<-relevel(ndata$dv2, ref="Vote")
over$ref<-relevel(over$dv2, ref="Vote")

# data for network victimization
no_direct_vic <- ndata %>%
  filter(victim_score_combined == "0")

no_direct_vic_over <- over %>%
  filter(victim_score_combined == "0")
#################################################################

#################################################################
# model 1 table (h1)
summary(model1)
pR2(model1)
n = nrow(fitted(model1))

stargazer(model1,
          title=c("Multinomial logit model (Model 1) testing the 
          relationship between victimization and agency attitudes. 
          Control variables for region are included in the model. 
          Voting is used as the reference category. Data: national sample."),
          covariate.labels=c("Low-level Vic.",
                             "High-level Vic.",
                             "East", "North", "Pacific", "South", "Intercept"),
          digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n)),
          out = "/Users/dorffc/Desktop/Replication/tables/model1_table.tex")
#################################################################

#################################################################
# model 2 table (h1 + controls)
summary(model2)
pR2(model2)
n = nrow(fitted(model2))

stargazer(model2,
          title = "Multinomial logit model (Model 2) testing the relationship between 
          victimization and agency attitudes. Control variables for age, gender, 
          education, familiarity with the Movement for Peace with Justice and Dignity, 
          and region are included in the model. Voting is used as the reference category. Data: national sample.",
          covariate.labels=c("Low-level Vic.",
                             "High-level Vic.",
                             "Age","Gender (woman)", "Education",
                             "Familiarity MPJD",
                             "East", "North", "Pacific", "South", "Intercept"),
          digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n)),
          out = "/Users/dorffc/Desktop/Replication/tables/model2_table.tex")
#################################################################

#################################################################
# model 3 table (h2)
summary(model3)
pR2(model3)
n = nrow(fitted(model3))

stargazer(model3,
          title = "Multinomial logit model (Model 3) testing the relationship between 
          network victimization and agency attitudes.
          Control variables for region are included in the model. 
          Voting is used as the reference category. Data: National sample subset.",
          covariate.labels=c("Network Victimization",
                             "East", "North", "Pacific", "South", "Intercept"),
          digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n)),
          out = "/Users/dorffc/Desktop/Replication/tables/model3_table.tex")
#################################################################

#################################################################
# model 4 table (h2 + controls)
summary(model4)
pR2(model4)
n = nrow(fitted(model4))

stargazer(model4,
          title = "Multinomial logit model (Model 4) testing the relationship between 
          network victimization and agency attitudes.
          Control variables for age, gender, 
          education, familiarity with the Movement for Peace with Justice and Dignity, 
          and region are included in the model. Voting is used as the reference category.Data: national sample subset.",
          covariate.labels=c("Network Victimization",
                             "Age","Gender (woman)", "Education",
                             "Familiarity MPJD", 
                             "East", "North", "Pacific", "South", "Intercept"),
          digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n)),
          out = "/Users/dorffc/Desktop/Replication/tables/model4_table.tex")
#################################################################

#################################################################
# model 5 oversample
# checks & set reference level
dim(over)
over$ref<-relevel(over$dv2, ref="Vote")

summary(model5)
pR2(model5)
n = nrow(fitted(model5))

# direct victimization
stargazer(model5,
          title = "Multinomial logit model (Model 5) testing the relationship between 
          victimization and agency attitudes using only data from the oversample. 
          Control variables for age, gender, and education level are included in the model. 
          Voting is used as the reference category.",
          covariate.labels=c("Low-level Vic.",
                             "High-level Vic.",
                             "Age","Gender (woman)", "Education",
                             "Familiarity MPJD", 
                             "Intercept"),
          digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n)),
          out = "/Users/dorffc/Desktop/Replication/tables/model5_table.tex")


#net victimization
summary(model6)
pR2(model6)
n = nrow(fitted(model6))

stargazer(model6,
          title = "Multinomial logit model (Model 6) testing the relationship between 
          network victimization and agency attitudes using only data from the oversample. 
          Due to sample size, controls are reduced.
          Voting is used as the reference category.",
          covariate.labels=c("Network Victimization", 
                             "Age","Gender (woman)",
                               "Intercept"),
          digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n)),
          out = "/Users/dorffc/Desktop/Replication/tables/model6_table.tex")
#################################################################

#################################################################
# model 7 table (h1 + alt. controls)
# robustness
summary(model7)
pR2(model7)
n = nrow(fitted(model7))

stargazer(model7,
          title="Multinomial logit model (Model 7) testing the relationship between
          victimization and agency attitudes including control variables age, gender,
          education, familiarity of the Movement for Peace with Justice and Dignity,
          urban/rural location, previous voting behavior, and region. Voting is used as the 
          reference category. Data: National sample.",
          covariate.labels=c(covariate.labels=c("Low-level Vic.",
                                                "High-level Vic.",
                                                "Age","Gender (woman)", "Education",
                                                "Familiarity MPJD", 
                                                "Rural", "Voted",
                                                "East", "North", "Pacific", "South", 
                                                "Intercept")), digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n )),
          out = "/Users/dorffc/Desktop/Replication/tables/model7_table.tex")
#################################################################

#################################################################
# model 8 table (h1 + alt. controls for SES)
# robustness
summary(model8)
pR2(model8)
n = nrow(fitted(model8))

stargazer(model8,
          title="Multinomial logit model (Model 8) testing the relationship between
          victimization and agency attitudes including control variables age, gender,
          socio-economic status, familiarity of the Movement for Peace with Justice and Dignity,
          and region. Voting is used as the 
          reference category. Data: National sample.",
          covariate.labels=c(covariate.labels=c("Low-level Vic.",
                                                "High-level Vic.",
                                                "Age","Gender (woman)", 
                                                "Num. of Lightbulbs",
                                                "Familiarity MPJD",
                                                "East", "North", "Pacific", "South", 
                                                "Intercept")), digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n )),
          out = "/Users/dorffc/Desktop/Replication/tables/model8_table.tex")
#################################################################

#################################################################
# model 9 table (political engagement)
summary(model9)
n = nrow(fitted(model9))

stargazer(model9,
          title="Multinomial logit model (Model 9) testing the relationship between
          victimization and agency attitudes with alternative controls for political engagement
          (control variables for party identification,
          information from community groups, and voting in the past election). Voting is used as the 
          reference category. Data: National sample.",
          covariate.labels=c(covariate.labels=c("Low-level Vic.",
                                                "High-level Vic.",
                                                "Party: PAN",
                                                "Party: PRD",
                                                "Party: PRI",
                                                "Info. Community Groups",
                                                "Voted",
                                                "East", "North", "Pacific", "South", 
                                                "Intercept")), digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n )),
          out = "/Users/dorffc/Desktop/Replication/tables/model9_table.tex")

#################################################################
summary(model10)
pR2(model10)
n = nrow(fitted(model10))

stargazer(model10,
          title="Multinomial logit model (Model 10) testing the relationship between
          victimization and agency attitudes with controls for participation in resistance groups
          (control variables for participating in protest and supporting or participation with
          self-defense groups). Voting is used as the reference category. Data: National Sample.",
          covariate.labels=c(covariate.labels=c("Low-level Vic.",
                                                "High-level Vic.",
                                                "Protest Participation",
                                                "Support/Participate in Self-Defense Groups",
                                                "East", "North", "Pacific", "South", 
                                                "Intercept")), digits=3,
          type = "latex", style="apsr", no.space=TRUE,
          add.lines = list(c("N", " ", " ", " ", n )),
          out = "/Users/dorffc/Desktop/Replication/tables/model10_table.tex")